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Abstract. We present calculations of quasiequilibrium sequences of irrotational binary neutron stars based on real- 
istic equations of state (EOS) for the whole neutron star interior. Three realistic nuclear EOSs of various softness 
and based on different microscopic models have been joined with a recent realistic EOS of the crust, giving in 
this way three different EOSs of neutron-star interior. Computations of quasiequilibrium sequences are performed 
within the Isenberg- Wilson-Mathews approximation to general relativity. For all evolutionary sequences, the in- 
nermost stable circular orbit (ISCO) is found to be given by mass-shedding limit (Roche lobe overflow). The EOS 
dependence on the last orbits is found to be quite important: for two 1.35 Mq neutron stars, the gravitational 
wave frequency at the ISCO (which marks the end of the inspiral phase) ranges from 800 Hz to 1200 Hz, depend- 
ing upon the EOS. Detailed comparisons with 3rd order post-Newtonian results for point- mass binaries reveals 
a very good agreement until hydrodynamical effects (dominated by high-order functions of frequency) become 
important, which occurs at a frequency ranging from 500 Hz to 1050 Hz, depending upon the EOS. 
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1. Introduction 

Detection of the gravitational waves (GW) with 
the use of the ground-based laser interferome- 
ters, such as V IRGO IjAcernese et al. 20040, L IGO 
H Abbot et al. 2004JI, GEO600 IjHewitson et al. 20030 . or 
TAMA300 l|Ando et al. 20010) will provide crucial infor- 
mation about various astrophysical objects. Among them, 
coalescing binary neutron stars (NS) seem particularly 
interesting in order to probe their interiors. Matched 
filtering techniques involving high-order post-Newtonian 
effects will provide the individual masses and spins of 
the NSs. The last few orbits before the final merger are 
dominated by the strong tidal forces acting between the 
components. The tidal deformation of NSs shape and of 
the matter distribution in the stellar interior are expected 
to depend rather strongly on the poorly known equation 
of state (EOS) of dense matter. The GW signal carries 
therefore some imprint of the EOS. In particular, the 
final frequency of the last stable circular orbit is strongly 
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correlated with the compactness parameter M/R 
l|Faber et al. 20021 |Taniguchi fc Gourgoulhon 2003| |, 

and thus can provide direct constraints on the 
theory of dense hadronic matter above the nu- 
clear saturation density. A related event potentially 
rich of informations about the dense matter EOS 
is the merger of neutron star and a black hole 
l|Prakash fc Lattimer 20041 IPrakash et al. 20040) . In 
addition to locating the transition from inspiral to the 
merger phase, computations of the last stable orbits of 
binary NSs are also required for providing initial data to 
compute the dynamical merger phase l|Shibata et al. 2 003 
and references therein). 

The last orbits of inspiraling binary neutron stars have 
been studied by a number of authors in the quasiequilib- 
rium approximation, and in the framework of Isenberg- 
Wilson-Mathews (IWM) approximation of general rela- 
tivity (see |Baumgarte fc Shapiro 20 03 for a review). The 
quasiequilibrium assumption approximates the evolution 
of the system by a sequence of exactly circular orbits (as 
the time evolution of the orbit is much larger than the 
orbital period). The IWM approximation amounts to us- 
ing a conformally flat spatial metric (the full spacetimc 
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metric remaining non conformally flat), which reduces the 
problem to solving only five of the ten Einstein equa- 
tions. Within these two approximations, the most real- 
istic studies have considered irrotational binaries, as op- 
posed to synchronized ones, for the viscosity of neutron 
star matter is far too low to ensure synchronization dur- 
ing the late stage of the inspiral (Bildst en fc Cutler~1 992 
IKochaneck 1992JI . 

The quasi-totality of the existing quasiequilibrium 
IWM studies (Bonazzolaet al. 1999, Maronetti et al. 1999, 
Uryu & Eriguchi 2000, Uryu et al. 2000, Gourgoulhon ct 
al. 2001, Taniguchi & Gourgoulhon 2002b, 2003) employ 
a polytropic EOS to model the neutron star interior. The 
only exception is the recent work of Occhsl in et al. 2004 
who have used two EOSs : (i) a pure nuclear matter EOS, 
based on a relativistic mean field model and (ii) a 'hybrid' 
EOS with a phase transition to quark matter at high den- 
sity. At 2 x 10 14 gem -3 (i.e. in the vicinity of nuclear 
density), these two EOSs are matched to a polytropic one 
with adiabatic index 7 = 2.86. This last assumption of 
Oechslin et al. is somewhat ad hoc, because the EOS of 
neutron star crust is very different from a polytrope, and 
its local adiabatic index is much smaller. Indeed, the crust 
polytropic index within the inner crust (which contains 
some 99.9% of the total crust mass) varies from 7 ~ 0.5 
near the neutron drip point to 7 ~ 1.6 in the bottom layers 
near the crust-core interface (Douchin & Haensel 2001). It 
is to be noticed that the crustal EOS plays a crucial role 
in defining the mass-shedding limit which marks the end 
point of quasiequilibrium binary configurations. 

In the present article, we study the last orbits of irrota- 
tional binary neutron star systems by using a set of three 
dense matter EOSs which are representative of the con- 
temporary many-body theory of dense matter. Contrary 
to lOechslin et al. 20041 we describe the neutron star crust 
by means of a realistic EOS obtained in the many-body 
calculations. As in the works mentioned above, we use the 
quasiequilibrium and IWM approximations. 

The paper is arranged in the following way: in Sect. 
[21 the differences between EOSs used in the computation 
are briefly described. Sect. is devoted to a description 
of the numerical methods used to obtain the quasiequilib- 
rium orbital sequences. In Sect.0]the results are presented, 
whereas Sect. \5\ contains discussion and final remarks. 

2. Description of the equations of state 

2.1. The EOS of the crust 

The outer layer of neutron star contains neutron rich nu- 
clei, which due to Coulomb repulsion form a crystal lattice 
if matter temperature is below the corresponding melt- 
ing temperature. This layer is called neutron star crust, 
and extends up to the density at the crust-core interface, 
p cc ~ 10 14 g cm~ 3 . The precise value of p cc is model 
dependent, and varies within (0.6 — 1.4) x 10 14 g cm~ 3 . 
The EOS of the crust is rather well established (for re- 
view, see, e.g., IHaensel 20 03 ) . In the present paper the 



EOS of the crust was composed of three segments. For 
densities smaller than 10 8 g cm~ 3 we used the EOS of 
|Baym et al. 1971| For 10 8 g cm -3 < p < 10 11 g cm -3 we 
applied the EOS of Hacns ~& Pichon 19941 which makes 
maximal use of the experimental masses of neutron rich 
nuclei. Finally, for densities 10 11 g cm -3 < p < p cc we 
used the EOS of IDouchin fc Haensel 200ll For neutron 
stars of M = 1.35 M , the crust contains less than 
2% of stellar mass. However, it is the region most eas- 
ily deformed under the action of the tidal forces resulting 
from the gravitational field produced by the companion 
star. Below melting temperature, elastic shear terms in 
the stress tensor are nonzero, but they are two orders of 
magnitude smaller than the main diagonal pressure term 
plaensel 2001(1 . Dissipative and thermal effects accom- 
panying matter flow inside neutron stars will be briefly 
discussed in Sect. |S] they are expected to be small at 
the quasiequilibrium evolution stage. Therefore, to a very 
good approximation, the crust layer behaves in the tidal 
forces field as an ideal degenerate fluid, described by a 
zero temperature EOS. 

An important quantity which actually determines the 
response of the crust layer to the compression or decom- 
pression of matter is the adiabatic index 7 = (n/p)dp/dn, 
where p is the local pressure and n the corresponding 
baryon number density. In the outer crust, the pressure 
is determined by the ultrarelativistic electron gas, so that 
7 = 4/3. However, the outer crust contains only 10 -5 
of the NS mass. Some thousand times more, i.e, about 
O.O3M0, is contained in the inner crust, composed of a 
lattice of heavy neutron rich nuclei immersed in an neu- 
tron gas and an electron gas. Its adiabatic index depends 
strongly on the density, varying from 7 ~ 0.5 near the 
neutron drip point (extreme softening of the EOS), up to 
7 ~ 1.6 close to p cc (Douchin & Haensel 2001). 



2.2. The EOS of the core 

Neutron star core contains matter of density significantly 
larger than the normal nuclear density, equal to the satu- 
ration density of infinite symmetric nuclear matter, po = 
2.8 x 10 g cm -3 , and corresponding to the baryon num- 
ber density no = 0.16 £m -3 . For p > po the EOS of the 
core is poorly known, and this uncertainty grows rapidly 
with increasing density. The theoretical EOSs, derived us- 
ing different theories of dense matter, and different meth- 
ods of solution of the many-body problem, differ signifi- 
cantly at 10 15 g cm -3 , characteristic of the central cores 
of neutron stars under study. In the present paper, we 
will limit ourselves to neutron star cores consisting of nu- 
cleons and hyperons, i.e., baryons whose properties are 
known from terrestrial experiments. The speculative case 
of exotic phases of dense matter (kaon and pion conden- 
sation, quark deconfinement), whose existence is - as for 
this writing - not substantiated neither by experiments 
nor by astronomical observations, will be considered in 
future publications. The most extreme case of hypothcti- 
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cal strange stars built of self-bound strange quark matter 
will also be studied separately. 

Three EOSs of the core were considered. Two of them 
may be considered as a soft and stiff extremes of the 
EOSs of matter composed of nucleons, electrons, and 
muons. The first of them, BPAL12, is of phenomenolog- 
ical type and can be considered as a soft extreme of the 
nucleonic EOS of NS matter (JPrakash et al. 19971 see also 
Bo mbaci 1995JI . The second one, APR, is based on precise 
variational calculations and includes realistic two-body 
(Argonne A18) and three-body (Urbana UIX) nucleon in- 
teractions IjAkmal et al. 19 98) 1 . We considered also one 
EOS in which hyperons are present at p > pHi where the 
threshold for the hyperon appearance pn ~ 2po (model 3 
of |Glendenning 19 85); it will be referred to as the GNH3 
EOS. This EOS was obtained using the Relativistic Mean 
Field model of baryonic matter. For p < pn (nucleons 
only) this EOS is very stiff but causal (w SO und < c). The 
appearance of hyperons strongly softens the EOS as com- 
pared to the pure nucleon case. The hyperons soften the 
EOS because they are more massive than nucleons, and 
when they start to fill their Fermi sea they are slow and 
replace the highest energy nucleons. 
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Fig. 1. The pressure p against the energy density p for 
the EOSs used in the paper: APR (dashed line), GNH3 
(dotted line), and BPAL12 (solid line). The empty cir- 
cles correspond to the central density of a non-rotating 
stellar model with a gravitational mass equal to 1.35 Mq 
(Tabled) 



All three EOSs are displayed in Fig. Q] As we men- 
tioned they are very different because they assume very 
different strong interactions models at supranuclear den- 
sities. 
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1 The APR EOS violates the causality at high density region 
which appears in NS core only close to the maximum mass of 
the neutron star. 
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Fig. 2. Gravitational mass of static isolated stars against 
areal radius for the APR (dashed line), GNH3 (dotted 
line), and BPAL12 (solid line) EOSs. 



2.3. Mass-radius relation for static NSs 

The differences among the three EOSs are reflected in the 
M(R) curves of static stable models shown in Fig. [3 (M 
being the gravitational mass and R the areal radius) . 

The BPAL12 EOS is overall quite soft, and there- 
fore produces NSs with R decreasing rapidly with in- 
creasing M. The BPAL12 EOS yields the smallest R(M) 
for M > 1.2Mq. The maximum allowable mass for this 
EOS is marginally consistent with observations, M max = 
1.46 M . The APR EOS is stiff, with radius staying nearly 
constant at R ~ 11.5 km for M = (0.6 — 1.6)M , and 
yields high M max = 2.2 M Q . The BPAL12 and APR EOSs 
are soft and stiff extremes within our set of the EOSs as 
far as the values of M max are concerned. However, their 
M(R) intersect at M = 1.2M Q . Therefore, these EOSs 
produce NSs with similar values of compaction parameter 
M/Riov M~ 1.2M . 

The GNH3 EOS is different from the two previous 
ones, because its M(R) curve is composed of two different 
segments. The lower-mass segment (M < 1.3 M Q ) con- 
sists of stars with no or with only small hyperon cores. 
The radius stays nearly constant at R ~ 14 km for 
M = (0.5 — 1.3)A/q. This segment is connected via a 
"'knee" with a high-mass segment consisting of neutron 
stars with increasingly larger soft hyperon cores. Along 
this high-mass (hyperon-softened) segment, R decreases 
rapidly with increasing M, reaching a very flat maximum 
at M~ 1.9 M Q . 

Different aspects of the EOSs show up if we consider 
the 1.35 M Q static stars. Their parameters, calculated us- 
ing three EOSs are given in Tabled Let us remind that a 
NS configuration is a functional of the EOS at the densi- 
ties p smaller than the central density p c . We see that for 
the GNH3 EOS central density is barely higher than the 
hyperon threshold pu, and therefore the EOS inside NS 
is very stiff, actually the stiffest of all three. The largest 
stiffness of the GNH3 EOS in the interior of a 1.35 M Q 
NS can be clearly recognized by looking at Fig. 1 and Fig. 



Bejger et al.: Nuclear EOS and the last orbits of BNS 



EOS 


M/R 


R [km] 


Mb[M @ ] 


p c [10 14 g/cm 3 ] 


GNH3 


0.140 


14.262 


1.45351 


6.26 


APR 


0.176 


11.350 


1.49110 


9.80 


BPAL12 


0.191 


10.447 


1.48472 


20.22 



Table 1. Properties of isolated static neutron stars of 
gravitational mass M = 1.35 M for the three EOSs used 
in our computations. M/R = GM/Rc 2 is the compaction 
parameter, R is the areal radius, Mb is the baryon mass, 
and p c is the central energy density, respectively. 



3. Looking at Fig. 1, we see that this EOS has the high- 
est P at any p < p c . Additional information is contained 
in Fig. 3: this EOS has particularly large 7 for p ~ po, 
i.e., in the outer layers of NS core, which are therefore 
quite "inflated" in comparison with those in the other two 
NS models. These features are responsible for particularly 
large R for the GNH3 EOS. The GNH3 configuration has 
definitely the largest R and therefore the smallest com- 
paction parameter M/R. The BPAL12 EOS is clearly the 
softest, with R smaller by 27% than for the GNH3 EOS. 
Final, the APR EOS, which is moderately stiff in this 
mass range, yields R which is only 8% larger than for the 
BPAL12 model. 

The differences in stiffness reflect the characteristics 
of the nuclear model underlying each of the EOSs of the 
NS core. Using the density dependent adiabatic index 7 
it is particularly easy to visualize these differences, see 
Fig. The strong drop in 7 above pn reflects the hyperon 
softening in the GNH3 EOS. The values of 7 increasing 
up to the maximum at p — p c tell us about the stiffening 
of the APR EOS at p ~ 2p , in contrast to the behavior of 
the GNH3 EOS which softens close to this density. Finally, 
the BPAL12 EOS remains close to a polytrope with 7 ~ 
2.2, except for a small region around pq. 

3. Numerical method 

3.1. Numerical code for close binary systems 

The present computations of close binary neutron star 
systems rely on the assumption of quasiequilibrium state 
(helical Killing vector approximation), with irrotational 
flow of the fluid, and conformally flat spatial 3-metric 
(Isenberg-Wilson-Mathews approximation). We construct 
the quasiequilibrium sequences of binary neutron stars 
described by the realistic EOSs using a numerical code 
which solves the five coupled, nonlinear, elliptic equa- 
tions for the gravitational field, supplemented by an el- 
liptic equation for the velocity potential of irrotational 
flows. The code has already been used successfully for 
calculating the final phase of inspiral of binary neu- 
tron stars described by the polytropic equation of state 
dTaniguchi et al. 200 1| |Taniguchi fc Gourgoulhon 2002a| 
2002b and 2003). This code is built upon the C++ 
library LORENE (http://www.lorene.obspm.fr) and 
can be downloaded freely from Lorene CVS reposi- 
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Fig. 3. The adiabatic index 7 against the energy density 
p for the EOSs used in our computations: APR (dashed 
line), GNH3 (dotted line), and BPAL12 (solid line). The 
empty circles correspond to the central density of a non- 
rotating stellar model with a gravitational mass equal to 
1.35 M Q . 



tory, as Lorene/Codes/Bin_star/coal.C. The complete 
description of the resulting general relativistic equa- 
tions, the whole algorithm, as well numerous tests of 
the code can be found in ( |Gourgoulhon et al. 200 1) . 
Additional tests have been presented in Sect. Ill of 
dTaniguchi fc Gourgoulhon 2003) . 

The numerical technique relies on a multi-domain 
spectral method with surface-fitted coordinates. We have 
used one domain for each star and 3 (resp. 4) domains 
for the space around them for a large (resp. small) sepa- 
ration. In each domain, the number of collocation points 
of the spectral method is chosen to be N r x Ng x N v = 
25 x 17 x 16, where N r , Ng, and N v denote the number 
of points in the radial, polar, and azimuthal directions 
respectively. The accuracy of the computed relativistic 
models was estimated using a relativistic generalization 
of the the Virial Theorem (Friedm an et al. 20021 see also 
Sec. III. A of ITaniguchi &: Gourgoulhon 2003) . The virial 
relative error was a few times 10 -5 . 



3.2. The velocity potential of irrotational flows 

Let us briefly discuss a technical difference between the 
computations of the realistic EOS irrotational binary NS 
and the polytropic case considered by Bonazzola et al. 
1999, Maronetti et al. 1999, Uryu & Eriguchi 2000, Uryu 
et al. 2000, Gourgoulhon et al. 2001, and Taniguchi & 
Gourgoulhon 2002b, 2003. The difference stems from the 
fact that the realistic EOS presented in Sec. El are given 
in tabulated form, and a certain thermodynamic coef- 
ficient C, (required in the computation of the velocity 
potential of the irrotational fluid flow, cf. Sect. II of 
|Gour goulho n et al. 2001| ) is not given explicitly by the 
tabulated EOS. 
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As already discussed in Sect. 2, we adopt the approx- 
imation of the perfect fluid for the form of the stress- 
energy tensor, and we represent the NS matter by a zero- 
temperature EOS. In view of this, we can use the Gibbs- 
Duhem identity: 



h 



P 



(1) 



where e = pc 2 is the proper energy density of the fluid, p 
is the pressure, h = (e+p)/(mhn) is the specific enthalpy, 
with n denoting the fluid baryon number density, and mh 
the mean baryon mass. 

Following Gourgoulho n et al. 2001| we write the equa- 
tion for the coefficient C, as 



c 



d(ln#) _ d(hxff) d(lnp) _ d(hxff) 

d(lnn) d(lnp) d(lnn) d(lnp) 



(2) 



where 7 is the adiabatic index. The quantity H := ln/i 
appears in the first integral of fluid motion and denotes the 
log-enthalpy of the fluid, and the value of d(lni/)/d(lnp) 
can easily be evaluated by means of Eq. Q which, as a 
strict thermodynamic relation, is exact. 

In practice, we used two methods to get the adiabatic 
index 7 = (n/p)dp/dn. In the first method we used an 
analytical formula for the adiabatic index being approx- 
imation of tabulated values of 7. In the second method 
the adiabatic index was obtained directly from a tabu- 
lated EOS by taking the derivative of the second order 
polynomial which goes through three consecutive (p, n) 
EOS points. The value of the dp/dn multiplied by the n/p 
ratio is evaluated in the middle point, and the resulting 
discrete values of 7 are ready to be interpolated, similar 
to the other quantities from the tabulated EOS. 

The second method was proved to be very robust, and 
acts as a consistency check; it was tested with EOS for 
which the precise values of the adiabatic index were ob- 
tained from the microscopic considerations, namely for the 
SLy EOS UDouchin fc Haensel 2000)1 . 

4. Numerical results 

4.1. Quasiequilibrium sequences of neutron stars 
described by realistic EOS 

In this section we present the numerical results for evolu- 
tionary sequences of close neutron star binaries described 
by the three realistic EOS introduced in Sect.El By evolu- 
tionary sequence, we mean a sequence of quasiequilibrium 
configurations of decreasing separation d and with con- 
stant baryon mass Mb- Such a sequence is expected to ap- 
proximate pretty well the true evolution of binary neutron 
stars, which is entirely driven by the reaction to gravita- 
tional radiation and hence occur at fixed baryon number 
and fluid circulation (zero in the irrotational case consid- 
ered here). We calculated evolutionary sequences of binary 
system composed of two identical neutron stars (equal 
mass system), with the baryon mass Mb corresponding 
to the gravitational mass M := Mi = M<x = 1.35M for 



a static isolated star (see Tabled for the values of Mb). 
We have selected the gravitational mass M = 1.35M Q 
for two reasons: (i) it agrees with "average NS mass" ob- 
tained from observations of radio pulsars in binary sys- 
tem and (ii) it allows us to compare our results with cal- 
culations made by other authors for polytropic models 
<|Faber et al. 2002)1 . 

For the benefit on the present discussion, let us re- 
call the definition of the ADM mass of the system, which 
measures the total energy content of a slice t — const of 
spacetime (hypersurface E 4 ): 



Madm '= 



1 
16tt 



[V^ id - V t (/ w 7w)] dS 



(3) 



where 7^ is the metric induced by the spacetime metric on 
St, fij is a fiat metric on £ t (the condition of asymptotic 
flatness being 7^ — > fij), T>i is the covariant derivative 
associated with /y, and the integral is to be taken on a 
sphere at spatial infinity. In the case considered here of a 
conformally flat spatial metric, 7^ = ty 4 fij, the surface 
integral J3J can be converted into a volume integral over 
the whole hypersurface £*: 



M A 



DM • — 



X, 



&\E+—K ii K l >\#x i 



(4) 



where E is the total energy density of matter with respect 
to the observer whose 4-velocity is normal to E t and Ky 
denotes the extrinsic curvature tensor of E t . 

At infinite separation, the ADM mass of the system, 
Eqs. J3J) and Q, tends toward the sum of the gravitational 
masses of isolated static stars, and will be denoted by M M : 



lim Madm = Moo ■=M 1 + M 2 = 2.7M©. 

d — >oo 

We then define the orbital binding energy by 
-Ebind := Madm — M^. 



(5) 



(6) 



The variation of -©bind along an evolutionary sequence cor- 
responds to the loss of energy via gravitational radiation. 
Gravitational waves are emitted mostly at twice the or- 
bital frequency: /gw = 2/. 

4.2. Innermost stable circular orbit 

Each evolutionary sequence terminates by a mass- 
shedding point, which marks the end of existence of 
quasiequilibrium configurations. The shape of the stars 
close to this limit is presented in Fig.0] The mass-shedding 
is revealed by the formation of a cusp at the stellar surface 
in the direction of the companion (Roche lobe overflow). 
This cusp is marginally visible in Fig. 0J 

A turning point of -Ebind along an evolution- 
ary sequence would locate an orbital instability 
(Fried man et al. 2002)1 . This instability originates both 
from relativistic effects (the well-known r — 6M last stable 
orbit of Schwarzschild metric) and hydrodynamical effects 
(for instance, such an instability exists for sufficiently stiff 
EOS in Newtonian regime, see e.g. |Taniguchi et al. 2001| 
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Fig. 4. Baryon number density isocontours in the coordinate plane y = (containing the rotation axis) (left panels) 
and in the coordinate plane z = (orbital plane) (right panels), for configurations close to the mass-shedding limit. 
The upper (resp. middle, lower) panels correspond to the GNH3 (resp. APR, BPAL12) EOS. The thick solid lines 
denote stellar surfaces. 



and references therein). It is secular for synchronized 
systems and dynamical for irrotational ones, as those 
considered here. Thus the quasiequilibrium inspiral of 
binary neutron stars can terminate by either the or- 
bital instability (turning point of -Ebind) or the mass- 
shedding limit. In the latter case, one may wonder 
whether there might exist equilibrium configurations be- 
yond the mass-shedding limit, i.e. dumb-bell like configu- 
rations (see e.g. |Eriguchi fc Hachisu 1985| ). However dy- 
namical calculations for 7 = 2 polytrope have shown 
that the time to coalescence was shorter than one or- 
bital period for configurations at the mass-shedding limit 
UShibata fe Urvu 2001|Marronetti et al. 2004|) . Therefore 
we may safely define the end of quasiequilibrium inspiral 
by the mass-shedding limit in the case where no turning 
point of -Ebind is found along the sequences (which is ac- 
tually the present case, as we shall discuss below). 

The variation of the orbital binding energy along evo- 
lutionary sequences is presented in Fig. [SI where the points 
correspond to the equilibrium configurations of binary sys- 
tem calculated using numerical method (see Sect. |3J) and 
the lines present our best fit described in detail in Sect. 4.3 
In the scale of Fig. there is no visible difference be- 



tween our numerical results and the fit, i.e. fitting curves 
pass through the points. We plotted also the binding en- 
ergy obtained in the framework of Newtonian theory and 
the 3PN post-Newtonian approximation for point masses 
( Blanchet 2002 ) . Figure shows clearly that no turning 
point of -Ebind occurs along evolutionary sequences. Hence 
there is no orbital instability prior to the mass-shedding 
limit. We conclude that the innermost stable circular orbit 
(ISCO) of the computed configurations are given by the 
mass-shedding limit. 

To compare different approaches to the relativistic de- 
scription of binary systems we present in Fig.Elthe orbital 
binding energy after subtraction of the Newtonian term 
oc Jq W - This figure shows thus the effects of general rela- 
tivity and of finite sizes (hydrodynamics). One can see that 
binary neutron star systems are quite far from Newtonian 
and 1PN systems. On the contrary, 2PN and 3PN results 
bv lBlanchet 20021 and 3PN results bv lDamour et al. 20001 
are very close to our results for a wide range of fre- 
quencies. Note that the difference between 2PN and 3PN 
approximations are much larger in IDamour et al. 2 000 
treatment (Effective One Body method) than in that of 
Blanchet 2002 (standard post-Newtonian expansion). 
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Fig. 5. Orbital binding energy -Ebind = Madm — M M 
of the binary system versus frequency of gravitational 
waves (twice the orbital frequency) along three irrota- 
tional quasiequilibrium sequences. The lines (dotted - 
GNH3, dashed - APR, solid - BPAL12) were plotted using 
the fit described in the text, whereas the points represent 
actual data. Thin solid line touching the bottom-right cor- 
ner shows the 3rd post-Newtonian approximation for point 
masses derived bv lBlanchet 20021 The lower thin dotted 
curve corresponds to the Newtonian limit for point masses. 
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Fig. 6. Orbital binding energy of the binary system mi- 
nus the (point mass) Newtonian term — fcjv/ 2 ^ 3 versus fre- 
quency of gravitational waves (twice the orbital frequency) 
along three irrotational quasiequilibrium sequences. Thin 
solid line shows the 3rd post-Newtonian approximation for 
point masses bv lBlanchet 2 002 Slightly below there is the 
2PN line (thin, dashed) and the 1PN one (thin, long dash- 
dotted). For comparison post-Newtonian approximation 
by IDamour et al. 20001 - 3PN (dashed-dotted thin line) 
and 2PN (thin, long dashed) are also plotted. 



4.3. Analytical fits to the numerical results 

Following IFaber et al. 20021 we have performed some 
polynomial fits to each of the computed sequences. This is 
a very important step in order to obtain the first deriva- 
tive of the functions required for the energy spectrum 
of gravitational waves (see below). We used two differ- 
ent approaches to approximate our numerical results. The 
first one, similar to that presented by IFaber et al. 20021 
is based on quadratic approximation of numerical results. 
First of all we decided to make a fit only to the differ- 
ence between obtained exact results for the ADM mass 
Madm [cf- Eq. (jHJ)] and the prediction of the Newtonian 

2/3 

theory - i.e. we made a fit to the function Sbind + ^n/gw' 
where the second term corresponds to the Newtonian 
point-like mass behaviour with &n = (G7r/4) 2 / 3 M 5 / 3 = 
4.06 10" 4 M Hz~ 2/3 . We found it sufficient to fit the nu- 
merical results by a second order polynomial without any 
linear term: 



-Ebind — -^n/qw 



fe/GW) 



(7) 



i.e. contrary to IFaber et al. 20021 we assume ki = 0. The 
best-fit coefficients h^ are collected in Table It should 
be mentioned that our aim was to obtain accurate for- 
mula for gravitational radiation in the region where it is 
effectively emitted and therefore we have performed the 
fitting procedure for frequencies larger than 500-600 Hz. 
This approximation works quite well for rather large fre- 
quencies (i.e. small distance between stars) and there is 
no need to introduce the linear term, as done by Faber et 
al (2002). The advantage of our quadratic formula is its 



simplicity and good accuracy in the region corresponding 
to the effective emission of gravitational waves. 

However it is possible to find much better approxima- 
tion of the numerical results if one takes into account high 
order post-Newtonian approximation for the binding en- 
ergy of point-mass systems. The 3PN formula as obtained 
bv lBlanchet 2 002 from the standard 2 post-Newtonian ex- 
pansion reads 



p3PN 
-^bind 

Moo 



1 L^ ~| i "* 



384 



1069 
3072 



il. 



3072 



41V - 



285473 \ 



864 



n 8/3 



(8) 



where 51* is the orbital angular frequency expressed in 
geometrized units: 

5Z* := Zn — / = 27T — rjew- 

The terms in 51» 2/3 , 51» 4/3 , 51* 2 and SI, 873 in Eq. © are 
respectively the Newtonian, 1PN, 2PN and 3PN term. 

In Fig. we present the difference between our numer- 
ical results and the 3PN approximation given by Eq. (JSJ. 
Formula (JSJ approximates very well the behavior of binary 
system of realistic neutron stars for a very large range of 
binary period (notice the scale of the y-axis of Fig. !). 
From Fig. we can define frequencies / np m as those fre- 
quencies at which the deviation from point-mass behavior 
becomes important. The values of these frequencies for 

2 i.e. non-resummed, in contrast to the so-called Effective 
One Body approach of Dam our et al. 2000 
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Fig. 7. Difference -Ebind — -^b™ between the binding en- 
ergy of irrotational binary neutron stars build upon real- 
istic equation of state and the binding energy of binary 
point masses in the 3PN approximation of Blanchct 2002 
The dots correspond to numerical results and the lines to 
polynomial fits to them (see text for details). 



each of the three considered EOSs are given in Tabic [3] 
Of course, they are not precise numbers and should be 
treated as having a ~ ±50 Hz uncertainty. We can ap- 
proximate the results presented in Fig.[7|by the power-law 
dependence on frequency /gw : 



Eh 



bind 



P 3PN _ „ 
^bind — a 



/< 



GW 



1000Hz 



(9) 



Because of the steep character of the function £kind— -^b™ 
seen in Fig. [7| the power n is quite large. The values are 
listed in Tabic It should be stressed that we assume 
the integer number of the power n, although this is in 
principle not obvious from the theoretical point of view. Of 
course the more careful treatment of this approximation is 
possible (eg. more terms in the expansion) but from Fig. 
it is clear that there is no need to do it and the leading 
term of high power is sufficient. 

One can draw an important conclusion from the pre- 
sented results and their comparison with relativistic ap- 
proximations for point masses in binary system. We can 
expect that taking into account next orders in post- 
Newtonian approximation does not change the energy by 
an amount larger than the difference between 2PN and 
3PN models. As a consequence the large deviation of our 
numerical results from the 3PN approximation is caused 
by the effect of a finite size of the star (e.g. tidal forces). 
The very high power n in relation © indicates that, even 
for small departures from point mass approximation, high- 
order tidal effects are very important, and dominate the 
relation -Et>md(/Gw)- Indeed the lowest order tidal term 
is known to be n = 4 (jLai et al. 1994(1 and the values ob- 
tained here are well above this. 



EOS 


k 2 [1O" 9 M s 2 ] 


a [AT©] 


n 


GNH3 


6.23 


2.745 KT a 


8 


APR 


5.38 


1.912 lO -41 


9 


BPAL12 


4.59 


6.519 KT 13 


16 



Table 2. Parameters of polynomial fits Q and @. 



4.4. Energy spectrum of gravitational waves for 
different realistic EOS 
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Fig. 8. Energy spectrum of GW waves emitted by the 
binary neutron star system versus frequency of gravita- 
tional waves (twice the orbital frequency) along three ir- 
rotational quasiequilibrium sequences. The straight lines 
correspond to the Newtonian dependence of energy mul- 
tiplied by 1, 0.8, 0.75 and 0.6. 

We computed the energy spectrum of gravitational 
waves obtained as the first derivatives of the fitted func- 
tions [Eq. ©]. The relation between d-Ebind/d/ and the 
GW frequency /gw is presented on Fig.|S| In this figure we 
draw also straight lines corresponding to the Newtonian 
case ~ /gw ^° nnc ^ the break frequencies at which the 
energy spectrum has dropped by 20%, 25%, 40%. These 
values are important from the point of view of the future 
detections: they show the difference between the ampli- 
tude of the real signal and the Newtonian template which 
allow to calculate the real wave form amplitude from the 
detector noise. We also compare our results with the 25% 
case from lFaber et al. 20021 It should be mentioned that 
there is no visible difference between our models for dif- 
ferent EOSs at the break frequency level of 10% (the case 
considered bv lFaber et al. 2002|) and the situation is then 
very precisely described by the 3PN formula. 

4.5. Comparison with polytropic EOS 

Up to now, all calculations (except those of 
IQechslin et al. 2004J) of the hydrodynamical inspiral 
and merger phases have been done for the simplified 
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EOS 


/npm 


ho 


ha 


/40 


J end 


GNH3 


500 


567 


657 


792 


806 


APR 


700 


615 


762 


1025 


1132 


BPAL12 


1050 


615 


785 


1160 


1270 



Table 3. Gravitational wave frequencies (in Hz) com- 
puted from the calculated data for GNH3, APR and 
BPAL12 EOSs: the / npm denotes the frequency at which 
the non-point-mass effects start to be important, /20, fm 
and /40 are the so-called break- frequencies (see text), 
whereas / en d is the GW frequency on the final orbit. 



equation of state of dense matter, for the poly tropic 
EOS, where the dependence between pressure and baryon 
density is given by p = nn 1 . It has been shown that 
the results obtained for these polytropic EOSs depend 
mainly on the compaction parameter M/R. It is therefore 
interesting to check if the properties of inspiraling neutron 
stars described by realistic EOS can be predicted, in a 
good approximation, by studying binaries with assumed 
polytropic EOSs. In order to make such a comparison we 
construct sequences of binary NSs described by polytropic 
EOS parametrized by compaction parameter M/R given 
in the Table Q] for each of three realistic equation of 
states. For a given 7, we calculate the value of the k 
coefficient which yields the same R at M = 1.35M Q as 
that predicted by a selected realistic EOS used in the 
present paper. The values of 7, the pressure coefficient k, 
and the baryon masses of polytropic static isolated NSs 
of M = 1.35M Q are presented in Table H 



corresponding EOS 


7 


k [pc 2 h 7 ] 


Mb [Mq] 


GNH3 


2. 


0.02645 


1.44802 


APR 


2. 


0.02037 


1.47016 




2.5 


0.01094 


1.49561 


BPAL12 


2. 


0.01908 


1.47742 



Table 4. Pressure coefficient k, adiabatic index 7 and 
baryon mass Mb for polytropic NSs having the same 
compactness parameter and mass (M = 1.35 Mq) than 
NSs described by realistic EOSs (Table QJ [p := 1.66 x 
10 14 g/cm 3 and h := 0.1 £m~ 3 ). 



In Fig. El we show the dependence of gravitational 
mass versus radius for isolated non-rotating stars based 
on the APR EOS, and two different polytropes having 
M/R=0.176 and M = 1.35M© at infinite separation. 
Finally in Fig. 1101 we present the binding energy versus 
the frequency of gravitational waves at the last orbits of 
the inspiral for the APR EOS and the two correspond- 
ing polytropes. The differences between quasiequilibrium 
sequences described by realistic and polytropic sequences 
are small. For the three different EOSs (APR and two 
polytropes) the frequency of gravitational waves at the 
last calculated orbit (/end) is ~ 1140Hz. Also the binding 
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Fig. 9. Gravitational mass versus stellar radius 
for sequences of static neutron stars described by 
lAkmal et al. 19981 EOS (solid line) and polytropic EOS 
with 7 = 2 (dotted line) and 7 = 2.5 (dashed line). 
Configurations with gravitational mass 1.35M Q (marked 
by a circle) described by polytropic EOSs have the same 
compaction parameter, M/i?=0.176, as a neutron star 
configuration based on the APR EOS. 
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Fig. 10. Binding energy of a binary neutron star systems 
as a function of gravitational waves frequency for the APR 
EOS (thick solid line) and two corresponding polytropic 
EOS with 7 = 2.5 (dashed line) and 7 = 2 (dotted line). 
Static isolated NSs have M/R = 0.176 and M = 1.35 M@. 
Thin solid line shows 3PN approximation for point masses 
bv lBlanchet"2002l Eq. ©. 



energy of the system at the mass-shedding point is close 
to each other, between —0.0372 and —0.0366 Mq. However 
one can see the differences in the frequencies / npm at which 
non-point-mass effects start to be important. The / np m 
has the smallest value for the polytrope with 7 = 2.5 and 
the highest for 7 = 2. Although the matter in the stel- 
lar interior is stiffer for the APR EOS (7 ~ 3, see Fig. 01, 
which is also clear from Fig.0 the APR curve lies between 
7 = 2.5 and 7 = 2 polytropes in FigED This is due to the 
relatively soft equation of state for the crust for realistic 
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NS models, which makes the response of the crust to tidal 
forces different from that of a polytrope with 7 = 2.5 or 
7-2. 

We obtained similar results by comparing two other 
sequences of realistic EOS - BPAL12 and GNH3 with cor- 
responding polytropic cases. As it turned out, the irrota- 
tional flow is weakly affected by the changes in the EOS 
of core, but it is expected that the differences should be 
seen in the merger phase. The outer layers of the star 
(those with subnuclear densities, i.e., the crust), which 
are properly treated in the present paper have influence 
on the properties of binary system at the last stages of 
inspiral. However the crucial parameter which determines 
the energy-frequency spectrum (energy per frequency) of 
emitted gravitational energy is M/R. 



5. Concluding remarks 

We have presented a set of evolutionary sequences of bi- 
nary neutron stars based on three selected realistic EOSs. 
These EOSs are based on modern many-body calcula- 
tions. Three baryonic EOSs of neutron-star core have been 
joined with a recent EOS of neutron-star crust, and in this 
way we obtained three different models of neutron star in- 
terior, from the surface to the stellar center. We restricted 
to models of neutron star core without exotic phases (me- 
son condensates, quark matter). In this way, the differ- 
ences between the core EOSs reflect the uncertainties in 
the existing theories of the interactions in nuclear matter. 

In the present paper we considered only those con- 
stituents of dense matter, which have been studied in labo- 
ratory. We did not consider here phase transitions to hypo- 
thetic exotic phases of dense matter, which were proposed 
by many authors, but which still remain hypothetical and 
speculative. Results obtained for the NS-NS binaries with 
exotic-phase neutron-star cores and realistic envelopes will 
be considered in our future publications. Similarly, the 
case of a binary involving strange quark stars built of a 
self-bound strange quark matter, will also be presented in 
a separate paper. 

We have computed quasiequilibrium sequences of ir- 
rotational NS-NS binary by keeping constant the baryon 
mass to a value which corresponds to individual gravita- 
tional masses of 1.35 M Q at infinite separation. For a long 
time of evolution of the binary system its binding energy is 
very accurately given by the 3PN post-Newtonian formula 
for point-mass system. However the departure from this 
3PN model at low binary periods has a quite abrupt char- 
acter, presumably due to high order tidal effects. The se- 
quences end at the onset of the mass transfer between the 
stars (i.e. when a cusp forms at the surface of the stars). 
This point defines the ISCO since no turning points of the 
binding energy has been found along the sequences, which 
would have revealed some orbital instability. The gravita- 
tional wave frequency at the ISCO is 800 Hz, 1130 Hz and 
1230 Hz, for the GNH3, APR and BPAL12 EOSs respec- 
tively. 



In a recent work based on the numerical integra- 
tion of the full set of time dependent Einstein equa- 
tions, IMarronetti et al. 20041 have located the dynamical 
ISCO by comparing the time evolution of quasiequilib- 
rium initial data at various separations (see also Fig. 14 
of |Shibata fc Uryu 200l| ). This defines the true ISCO, as 
opposed to the "quasiequilibrium" ISCO obtained here. 
For a polytropic EOS with 7 = 2 and a compactness 
parameter M/R = 0.14, they obtain the ISCO at an 
orbital frequency which is 15 % lower than the mass- 
shedding frequency of the quasiequilibrium sequence of 
Taniguchi & Gourgoulhon 2002b. This makes us confident 
that the values of the gravitational wave frequencies given 
above are quite close to those of the end of the true inspi- 
ral. 

When comparing our results with those of the recent 
work of Oe chslin et al. 20041 one should stress the basic 
difference in the EOS of matter at subnuclear densities. 
Oechslin et al. represented the EOS of the crust by a poly- 
trope with the adiabatic index 7 = 2.86. In this way, they 
made their EOS quasi-continuous at the crust-core inter- 
face, because their nuclear core EOS is very stiff. However, 
such an EOS for the crust is very unrealistic, because the 
real 7 can be as low as 0.5 near the neutron drip point. 
As the crust EOS is crucial for the last stable orbits of the 
NS binary, this implies differences between theirs and our 
results, even if for a M = 1.35M model, the compactness 
M/R resulting from the GNH3 EOS and the ir EOS are 
very close. In particular, IQechslin et al. 20041 have found 
a turning point in the binding energy along their se- 
quences, resulting in a quasiequilibrium ISCO. This dif- 
ference is certainly due to the stiffness of the polytropic 
EOS used by them in the outer layers of the star. Let 
us recall that for polytropic irrotational binaries a turn- 
ing point ISCO exists only if 7 > 2.5 ( |Uryu et aT~2 000 
ITaniguchi fc Gourgoulhon 2003| ). 

In our calculations we treated the neutron-star matter 
as an ideal fluid. In other words, we neglected the elastic 
shear (in the crust - if not molten) and viscous (in the crust 
and in the core) terms in the matter stress tensor. There 
terms are believed to be small, but they lead to specific 
physical phenomena. In particular, the matter flow in NS 
interior will break the beta equilibrium between baryons 
and leptons, and this will imply a neutrino burst at the last 
stage of inspiral (see Haensel 1992). Moreover, dissipative 
processes will heat the matter. Both effects will be studied 
in future publications. 

All numerical results presented here, includ- 
ing the full binary configurations, are available 
to download from the LORENE main server 
http : //www . lorene . obspm . f r/data/ 
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